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Abstract 

The surface composition of charged Lennard-Jones clusters A^ + , composed of N particles (55 < N < 
1169) among which n are positively charged with charge q, thus having a net total charge Q = nq, is in- 
vestigated by Monte Carlo with Parallel Tempering (MCPT) simulations. At finite temperature, the surface 
sites of these charged clusters are found to be preferentially occupied by charged particles carrying large 
charges, due to Coulombic repulsions, but the full occupancy of surface sites is rarely achieved for clusters 
below the stability limit defined in this work. Large clusters (N = 1169) follow the same trends, with a 
smaller propensity for positive particles to occupy the cluster surface at non-zero temperature. We show 
that these charged clusters rather behave as electrical spherical conductors for the smaller sizes (N < 147) 
but as spheres uniformly charged in their volume for the larger sizes (N = 1169). 
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I. INTRODUCTION 



The question of ion localization at interfaces has aroused considerable interest in the past 
decade, from the study of salts, acids, and bases at liquid interfaces 1 , to nanodroplets 1 , atmo- 
spheric aerosols 2 , and hydrophobic droplets 3 . The affinity of ions for the very surface of bulk 
water is a subtle topic, where the specific properties of ions and their effect on the structure of the 
solvent come into play^-. However, in finite-sized charged clusters or droplets, the distribution 
of charges is further affected by the fact that they are confined to a roughly spherical volume. In 
such systems, one issue is the competition between charge localization inside or at the surface. It 
is this latter effect that we study in the present work, employing a model whose liquid structure is 
far simpler than that of water. One of the most famous theoretical problem related to this question 
is the determination of the configuration that minimizes the Coulombic energy of N point charges 
located on a sphere, as initially proposed by J. J. Thomson. Such a fundamental problem finds 
applicability in modern investigations such as the formation pattern of virus capsids 7 , fullerene 
stability 8 , or multielectron bubbles at helium interfaces 9 , and has been treated by a number of the- 
oretical approaches among which simulated annealing 10 , genetic algorithms 11 , basin-hopping'i 2 ^, 
and constrained global optimization 14 . 

Electrospray Ionisation (ESI) in Mass Spectrometry— is a domain where charged droplets are 
formed and some interrogations on the charges localization still persist. In these experiments, 
charged droplets composed of solvent and analyte molecules (typically biomolecules) are pro- 
duced at finite temperature. After several cycles of evaporation, the analyte carrying part of the 
total droplet charge may eventually be ejected from the droplet and be analysed by mass spectrom- 
etry techniques in the gas phase. Beside the open question on the fragmentation mechanisms of 
these charged droplets^r— , one important issue is the localization of the embedded charges, i.e. 
within or at the surface of the droplets, and the influence of this charge localization on the droplet 
stability. Recent atomistic simulations have shown that the common belief that positive charges 
are located at the droplet surface should be taken with caution^ .. Localization of the embedded 
charges within charged droplets and the subtle effects on charges derealization are the subjects 
tackled in the present investigation. 

Inspired by the work of Jortner and co-workers, we have modeled ESI droplets as charged 
clusters 21 with a simple and generic representation of the interactions between individual parti- 
cles. In the rest of this work, droplets containing a positive net total charge will be considered, 



2 



thus corresponding to positive-ESI in experiments. Conclusions will be transferable to negative- 
ESI mode. In the positive-ESI mode, droplets with a net total positive charge are mainly produced 
because of the voltage applied in the capillary interface that drives positively charged entities to- 
wards the vacuum chamber. This does not preclude that negative charged particles are present 
within the droplet, up to the point that pure water droplets electrosprayed in positive-ESI mode 
can form negatively charged droplets, and the reverse, i.e. pure water droplets electrosprayed in 
negative-ESI mode provide a majority (>70%) of positively charged droplets, as shown by Jarrold 
and coworkers^—. Note also that the droplets are electrosprayed from a mixture of solvents (wa- 
ter, methanol and their mixture are typically employed in ESI of biomolecules), biomolecules and 
counterions (salts and buffers from the initial solution), for typical applications in proteomics and 
genomics. Residual counterions are thus present in the final positively charged droplets formed. 
Other sources for the presence of negative charges within a positively charged electrosprayed 
droplet are typically collisions with metal surfaces (from the walls of the capillary interface and 
skimmers for instance) or electric fields above the Taylor limit^. These remarks convey the idea 
that positively charged droplets of interest to our present modeling will be composed of charged 
particles carrying a positive charge q > 0, but it is not unreasonable to model charged droplets 
with a total net positive charge carrying both positive q > and negative q < individual charges. 

In the present investigation, we model charged droplets as clusters composed of N indi- 
vidual particles among which n are charged particles, each one carrying the charge q (q > in 
most calculations, but we will also investigate charged clusters composed of both q > and q < 
charges). The cluster has a total net positive charge Q = nq. We investigate the influence of 
structural changes (ie, cluster size N, number of embedded charged particles n) as well as electro- 
static (ie, charge sign and value, polarisability) and thermodynamic (ie, temperature) changes on 
the propensity for the charged particles to occupy the surface of the charged clusters. We are 
interested in clusters composed of 55 < iV < 1169 particles. 

The present investigation is based on Monte Carlo with Parallel Tempering (MCPT) simula- 
tions. We use a simplified dimensionless form of the interaction potential between the particles, 
that could be qualified in modern language as a mesoscopic highly-coarse-grained representation. 
It is based on Lennard- Jones and Coulomb pairwise interactions between individual spherical par- 
ticles, completed by the addition of iV-body polarisation. The interactions are written in dimen- 
sionless units. Each individual coarse-grained particle is either neutral (q = 0) or charged (q 0), 
and may carry an induced dipole moment if polarisation effects are taken into account in the cal- 



3 



culation. Within this model, one coarse-grained particle can be seen as either one biomolecule (for 
instance one amino acid, one peptide, or even one protein, for the proteic family), or as one coarse- 
grained solvent particle mimicking a few microscopic solvent molecules. In this latter case, our 
representation follows the coarse-grained models of water developed in the recent literature, where 
one chargeless coarse grain can be mapped to up to four individual microscopic water molecules. 
See for instance the review by Klein and coll. 26 and the popular Martini coarse-grained model—. 
The neutral spherical particles in our model can thus be seen as coarse grains for the solvent, i.e. 
water or methanol (and mixtures of) as typically used in ESI of proteins and non-covalent com- 
plexes of protein o 2Si29 . The charged coarse grains mimic biomolecules, and there is no explicit 
charge transfers between solvent grains and biomolecules grains in the model (i.e. the ESI charg- 
ing process is not directly modeled). It is easy to think of such spherical coarse grains in terms of 
charged amino acids (either for the native charged residues or after protonation of neutral residues 
in positive-ESI) or in terms of charged folded peptides and globular proteins. Similar spherical 
coarse grains are also used in the literature in order to simulate highly coarse grained non-globular 
proteins, i.e. patchy proteins 30-32 where the patches account for the anisotropy of interactions (and 
thus non spherical shapes). Here we consider a simple isotropic coarse-grained model of proteins, 
but one improvement will certainly be to introduce such patches. Obviously, the charged grains can 
also represent ionized water (a grain of water including one H 3 + molecule in positive-ESI), or a 
residual counterion (a grain would be seen as the counterion and its immediate solvation sphere). 
Our purpose is not to develop coarse-grained models dedicated to a specific solvent or to specific 
biomolecules. Rather we want to use a generic representation of the interactions between individ- 
ual grains, and to that end dimensionless units are used. This is an advantageous representation as 
one single simulation can be used for describing several different cases, i.e. encompassing several 
solvent and analyte compositions, provided an inverse mapping between dimensionless units and 
real units is eventually done. 

Section|n]of this paper describes the dimensionless potential energy surface used in the present 
work and computational details about the Monte Carlo with Parallel Tempering method. Results on 
charge localization and its relation to the excess charge de localization are reported in Section [III] 
Concluding remarks are collected in Section HVl 
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II. METHOD 



A. Interaction model 



Our A^ + clusters (i.e. composed of A" spherical particles, n are positively charged with the 
individual charge q, total net charge of the cluster Q = nq) are bound by Lennard- Jones and 
electrostatic interactions (charge-charge, charge-induced dipole, induced dipole-induced dipole) 
between individual spherical particles. Polarisation is expressed according to the point dipole 
model of Thole 3 ^. Defining the usual Lennard- Jones reduced energy V* = V/u and length 
r* = r/a, where u and a are the Lennard- Jones pair well-depth and diameter parameters, the 
reduced charge and polarizability are respectively q* = qj '(47re a w) 1 / 2 and a* = a/(47re a 3 ), 
and the reduced dipole moment is p* = p/ (Aire^a 3 w) 1//2 , where eo is the vacuum permittivity. The 
dimensionless potential used in the present simulations is: 

= E 4 (^2-^) + s£tf Ml" 1 *? + ~E^+£*V;)0z* -P*Vmr* 3 ) (1) 

considering that all spherical particles (or grains) have the same size. In this equation, r*- is the 
dimensionless distance between particles i and j, q*, a* and p* are respectively the dimensionless 
charge, polarisability and induced dipole of particle i. 4> s (r*j) is the Thole damping function 



(40 = k 



i _ ( i + ^Jii \ e -«7« 



(2) 



with a = 2.1304 and 7^ = r* ? /(a*a*) 1 / 6 . The dimensionless polarisability is chosen equal to 
a* = 0.1 for all the polarisable particles and the charge q* = q* is the same for all the charged 
particles, unless specified otherwise (ie, when we consider droplets with positive and negative 
charged particles). 



B. Monte Carlo with Parallel Tempering (MCPT) simulations 

The MCPT simulations were performed in spherical containers of radii R = 4.5a, 5.0a, 5.5a, 
and 10a respectively for the cluster sizes A^ = 55, 100, 147, and 1169 (note that these values 
exceed the typical radii used for modelling the corresponding neutral clusters so as to provide 
more flexibility for geometrical distortion of charged clusters). The centre of mass of the cluster 
and the container centre are constrained to coincide. 
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Each parallel MCPT simulation is a set of N rep Monte Carlo simulations performed at different 
dimensionless temperatures T* with different initial random seeds £j. The temperatures T* encom- 
pass the temperature range from the solid-state (low-temperature asymptote of the heat capacity 
curves) to liquid-state temperatures (high-temperature asymptote of the heat capacity curves). We 
adjusted the temperature scale for each investigated cluster in order to get the better compromise 
between resolution of the melting peak and convergence of the simulation, maintaining the re- 
quired overlap between adjacent energy distributions between 40% and 90%. 35 

Three MC moves have been adopted in the MCPT simulations: (i) displacement of particles, (ii) 
swap between a neutral and a charged particle, (iii) swap between configurations from adjacent 
replicas at temperatures T* and T* ±1 . Simulations were performed for 10 6 to 10 8 sweeps (one 
sweep = iV MC steps) depending on the cluster size and whether polarisation is taken into account 
(more time consuming calculations). This number of MC steps and thus the convergence of the 
MCPT calculations have been checked on the heat capacity curves of neutral clusters, comparing 
our results to available data from the literature.- 3 -^ Global geometry optimisations have also been 
performed with the GMIN code 38 . 

C. Stability criterion 

Despite the cohesion due to Lennard-Jones interactions and polarisation, charged clusters 
should become rapidly unstable as the total reduced charge Q* = nq* increases. We therefore 
defined a stability limit below which clusters can be regarded as stable and are included and con- 
sidered in our analyses. The stability criterion is established at the melting temperature, and our 
definition is based on the fraction of configurations in which evaporation is detected. An individ- 
ual configuration is considered to undergo evaporation if at least one particle reaches the container 
edge within a small distance of 1/2 in reduced units (or a/2 without reduced units). For N < 147, 
the neutral clusters are not prone to evaporation at the melting temperature and a charged cluster 
ffme/i rou ghly in the range 0.23-0.36 for the charged clusters considered here) was considered 
stable if no evaporation occurred at this temperature among at least 90% of the configurations 
explored by the MCPT simulation. For N = 1169, the neutral cluster is already prone to some 
evaporation (in ~ 20 — 25% of the configurations) at the melting temperature and a charged clus- 
ter (T^ ielt roughly in the range 0.52-0.55 for the charged clusters considered here) was considered 
stable if no evaporation occurred at this temperature among at least 60% of the configurations 
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explored by the MCPT simulation. 

HI. RESULTS 

A. Localization of charged particles and temperature effects 

At T* = 0, the charged particles of A^ + clusters are located at the cluster surface. For in- 
stance, the global minimum configurations of icosahedral A\^ + and clusters have their twelve 
charged particles located on the twelve icosahedron vertices whatever the value of the charge q*. 
However, when the temperature increases, the global minimum configuration is (not unexpectedly) 
not systematically entropically favoured. This can be seen in Figured] where the energy difference 
AV* between the local minima belonging to the basins explored during the MCPT simulations 
and the global minimum is plotted for selected small clusters at solid-state (cyan), melting (black), 
and liquid-state (red) temperatures, respectively. Note that the melting temperature is defined as 
the temperature at the heat capacity peak, the solid-state and liquid-state temperatures are the min- 
imum and maximum temperatures of our temperature ranges respectively, and each distribution 
was obtained by performing local minimisations with a conjugate gradient method over at least 
10 3 configurations collected along the MCPT simulations. 

At low temperature, the energy distributions of the smaller clusters such as Al® + are located 
in the immediate surrounding of the global minimum (AV* = 0, FigJT]), all the more when q* 
is large. Conformations structurally close to the global minimum are explored and the charged 
particles are thus preferentially located at the cluster's surface, though possibly occupying sites 
on the faces or edges rather than the icosahedron vertices. For N = 147 with n = 74 embedded 
charges, the surroundings of the global minimum can only be explored for small q* values. Going 
from g* = 0.1 to q* = 0.5, AV* already amounts to 5-20, and the global minimum energy basin is 
no longer significantly explored during the MCPT simulation for q* = 0.5. Fully charged clusters 
such as A\^ + are peculiar cases with much simpler energy landscapes, preserving the icosahedron 
structure of the neutral cluster, and therefore exploring very few minima other than the global one. 
At liquid-state temperatures (red plots in FigfT}, the global minimum surroundings are no longer 
explored, and much higher energy conformations are found for all the systems presented in Fig. \T\ 
At these temperatures, the clusters adopt a variety of distorted geometries in order to maintain 
their ability to accomodate, to some extent, charged particles at their surface as discussed below. 
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We have plotted in Figures |2^-d the average fraction F sur f of charged particles located at the 
cluster surface as a function of the total cluster charge Q* (Q* = nq* for A^ + , Q* = (n — m)q* 
for A^ +m ~), at the melting temperature of the charged clusters. The fraction F surf is defined by 
the ratio of the number of charged particles located at the surface over the maximum number of 
charged surface particles (the method for surface particle definition is described in Ref. 39 ). A value 
of F surf = 1 would correspond to charged particles only located at the cluster surface. Note that 
the number of surface sites for neutral A 55 , Ai o, Ai 47 , An 69 clusters are respectively 42, 69, 92, 
and 449 at T* = 0. These values are subject to a 5-10% increase at the highest temperatures of the 
MCPT simulations. 

We first discuss the results obtained without polarisation for the smaller clusters, ie iV < 147 
(black open symbols in Fig.[2^-c). For these systems, F surf varies between 0.66 (see Fig. |2];) and 
0.99 (see Fig. [2^), at the melting temperatures. The propensity of the n positive particles to oc- 
cupy the surface of A^ + clusters steadily increases as Q* increases. This is due to the increase 
in the repulsive Coulomb interactions that repel the charges towards the cluster surface as q in- 
creases. F sur f sensitively depends on the repartition (n, q*) for a given cluster size N. Two clusters 
with equal Q* such as A™ + (q* = 0.7) and Af 5 + (q* = 0.5) [Q* = 21], or A}°+(<f = 1.0) and 
Al°+(q* = 0.5) [Q* = 10], or A{${q* = 1.0) and A^+( g * = 0.3) [Q* « 14], present different 
fractions of charged particles at the cluster surface: for a given Q*, stable clusters with a larger 
number of charged particles n have smaller charge values q* and therefore display smaller frac- 
tions of charged particles occupying the cluster surface sites. The upward error bars in Figs |2k-b 
correspond to F sur f values spanned at the lower temperatures of the MCPT simulations (solid- 
state). It shows that F sur f increases by 5 — 10% at solid-state temperatures, possibly reaching unity 
for certain cases of N = 55 and iV = 100. Temperature is thus a pivotal parameter on which the 
fraction of charged particles located at the cluster surface strongly depends. 

Figure |2] also reports F surf values when including polarisation, showing a slight decrease in 
the number of charged surface particles. For instance, F sur f decreases from 0.974 to 0.906 for 
AH + (q* = 1), or from 0.880 to 0.845 for A^(q* = 0.5). This loss of charged surface particles 
is significantly reduced for mixed clusters (ie, composed of both positive and negative charged 
particles). For instance, F surf only decreases from 0.716 to 0.710 for A^ 14 ~(g* = 0.5). 

For N = 1169, the general trends observed on small clusters, i.e. increase of F surf with Q* 
at fixed n and faster increase of F surf for small n at fixed Q*, are preserved. However, unlike 
smaller clusters, a smaller fraction of the charged particles are located at the surface of the A™^ 
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clusters. The ratio of surface sites compared to inner sites is lowered in larger clusters, and the 
full occupancy of surface sites is never achieved although n may greatly exceed the number of 
accessible surface sites. For instance, the number of surface sites of A™69 clusters is about 450 — 
490 over the temperature range spanned in the MCPT simulations but the cluster surface is never 
fully occupied (~ 72% maximum occupancy found for the cases investigated here). 

It is also relevant to explore whether this surface propensity for the positively charged parti- 
cles remains when both positive and negative charged particles are embedded within the cluster. 
We have investigated that point in the case of the prototype A 44 ^ 30 ~ and A*^ 14- clusters, both 
composed of a total of 74 embedded charged particles with fixed |g*|, respectively with 44 q* > 
charges and 30 q* < charges (A 44 | 30 ~, Q* = + 14|<f|) and 60 q* > and 14 q* < charges 
(A^ 14 ~, Q* = +46 1 5* |). These mixed clusters have therefore the same total charge as A} 4 ^ and 
A 4 47 respectively, and the same number of charged particles as A 7 ^ , but their fraction of charged 
particles on the cluster surface is significantly smaller (see Fig. |2fc, filled black symbols). The 
negative particles are more prone to reside inside the cluster than the positive particles. For in- 
stance, only one third of the negative particles present in Af°7 14 ~(|g*| = 0.5) clusters are located 
on the surface. Due to attractive Coulombic interactions, the positive particles tend to organise 
themselves around the negative particles, thus decreasing their presence at the clusters surfaces. 
There is thus a competition between solvating the negative particles within the cluster and reaching 
the cluster surface. As an example, we found that F surf w 0.81 for A^y 14 ~(|g*| = 0.5) which is 
below the fractions obtained for A 4 ^ (q* = 0.5) and A 7 ^ (q* = 0.5), respectively 0.85 and 0.88. 
Including polarisation does not change these conclusions. 

B. Are charged droplets spherical conductors or uniformly charged spheres? 

In a recent paper, Ahadi and Konermann 19 very nicely showed that although excess charges 
(Na + ions in their case) embedded in a water droplet are not located at the droplet surface the 
excess charge is confined to a thin layer at the droplet periphery, so that the droplet acts as an ideal 
spherical conductor. Their demonstration is based on the electrostatic potential energy between 
the charged droplet and an external point charge probe. The structures of the droplets are extracted 
from room temperature microscopic molecular dynamics simulations of 10 Na + ions embedded in 
a water droplet of roughly 20 A (~1250 water molecules), close to the Rayleigh limit. They found 
that the dependence of the electrostatic potential energy with respect to the distance r between the 
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probe and the center of mass of the droplet is identical to the one obtained for a spherical conductor, 
ie, constant within the radius of the droplet and varying as 1/r outside the droplet. The argument 
used in the paper in order to reconcile the apparent contradiction between the charges location 
(inside the droplet) and the excess effective charge confined to a thin layer at the droplet surface 
is based on charge-induced orientation of the water dipole moments: the water molecules solvate 
the Na + ions within the droplet which induces an "orientational polarization" or "orientational 
alignment" of the water molecules. This consequently acts in a way so that the excess charge 
is transfered to the droplet periphery. The fact that the water molecules are polar molecules that 
form intermolecular hydrogen bonds is probably pivotal in the final charge-induced orientation of 
the water dipole moments. One might then logically wonder how much of these arguments are 
transferable and applicable to droplets not composed of explicit microscopic water molecules, but 
instead composed of solvent molecules having neither a permanent dipole moment neither forming 
hydrogen bonds. Our model for coarse-grained charged droplets provides an excellent test, and we 
apply now the same electrostatic calculations as in Ref.— on selected charged clusters investigated 
in Section IIII Al Note that some calculations presented below will include polarization effects in 
terms of induced dipoles on the particles, which will account for "charge transfers" between the 
particles of the clusters (an effect not taken into account in Ref.—). 

Following Ahadi and Konermann 19 , we present in Figj3]the Coulomb potential energy U\ lc 
curves obtained from our MCPT simulations: 



between a point charge probe q* be = q* and all A" particles q* = q* of the charged cluster 
A^ + . In this equation, r* and r * are respectively the positions of particle i and of the probe point 
charge defined with respect to the cluster center of mass. Ul IC (r*) is averaged over at least 10 3 
Monte Carlo configurations, and anticipating that A^ + clusters are not exactly spherical we have 
furthermore averaged over the three directions of space and over three planes (xy, yz, zx). In 
Ref. 19 the droplet charge was close to the Rayleigh limit, we have therefore used some of our 
results obtained for clusters with a relatively high total charge Q*, although this charge may still 
be far from the Rayleigh limit. 

The electrostatic arguments behind Eq.© are now the following. If we consider our charged 
clusters as roughly spherical in shape, they might either behave as a spherical conductor charged 
at the surface or as a sphere uniformly charged in volume. Following the Gauss law for the 




(3) 
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calculation of the radial electrostatic field E*(r*) and the related electrostatic potential V*{r*) 
created by such objects anywhere in space, i.e. for r* < Rq and r* > Rq where Rq is the radius of 
the sphere, we get the following expressions for the Coulomb potential energy (see Appendix for 
the derivation of the expressions in reduced units): 

r* <R* : U* surf = 47iq* probe a* Q R* (4a) 
r* > Rl : U; urf = 4 ^;W 2 . (4b) 

for the spherical conductor charged at its surface with the surface density of charge ctq = Q* /S*, 
and 



^QprobePQ 



r 



*2 



r*<R* : U; ol = ^ * (iff - — j (5a) 

r* > : U; ol = A ^^ QRf . (5b) 

for the sphere uniformly charged in volume with the volume density of charge p*q = Q* /V*. The 
Coulomb energy between a spherical conductor or a sphere uniformly charged in volume and an 
external probe particle varies as 1/r* for r* > Rq in both cases. However, the behaviour largely 
differs for r* < Rq since this energy is constant (the electrostatic potential is constant) inside the 
spherical conductor and decreases as — r* 2 when the sphere is uniformly charged in volume. 

We have plotted in Fig. [3] the Coulomb energies U* MC extracted from the MCPT simulations 
(solid curves) at solid-state (cyan curves) and liquid-state (red curves) temperatures. We compare 
these results to the ideal Coulomb energies (dashed curves, same convention of colours as above 
for solid/liquid states) U* urf defined by Eqs. (l4al)-(l4bT) and U* ol defined by Eqs. (I5ab-(I5B1). Rq val- 
ues have been fitted so as to provide the best agreement with the curves from the MC simulations 
in terms of absolute Coulomb energy values. Note that the radii thus obtained agree very well with 
the ones anticipated from radial distribution functions of the clusters (not shown here). We might 
anticipate some deviations in the absolute values of U* MC and U* ur j (respectively U* ol ) since the 
clusters investigated here are not perfectly spherical, although some of them might have nearly 
icosahedral symmetry, but we will see below that deviations (when they exist) are very small in 
most of the cases investigated here. 

For iV < 147, the fraction of charged particles located at the cluster surface always exceeds 
~ 0.7 (see Fig. [2]), so that a majority of the charged particles are located at the cluster surface. We 
find that these clusters behave as ideal electrical spherical conductors at solid-state temperatures 
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where we can see that the Coulomb energy V* MC closely follows the curves calculated for a sphere 
uniformly charged at its surface (dashed cyan curves). The ideal curves were obtained with Rq = 
1.98, 2.50, and 2.69 for clusters with N = 55, 100, and 147, respectively (see Figs [3h-c). We 
observe deviations from the ideal behaviour for A^J" (q* = 0.5, Q* = 37) when r* = 2 — 3 
(close to the surface), which we assume are related to distortions from sphericity for these clusters 
because of the higher value of the total charge Q*. At liquid-state temperatures (red lines in Fig. |3), 
the clusters display a behaviour between spherical conductors and spheres uniformly charged in 
volume. The Coulomb energy curve is close to the spherical conductor case for r* < 2 where the 
potential energy is roughly constant (as expected for a conductor). Note that the constant is either 
equal to the one obtained at the solid-state temperature (A^qq ) or slightly lower than this reference 
value (Af 5 + and A^+). For r* > 2, the curves then depart from the ideal constant of the spherical 
conductor, but they nonetheless do not coincide with the ideal law of a sphere uniformly charged 
in its volume (dashed red curves which always stand below the red solid curves). The deviations 
from the ideal spherical conductor might therefore again be related to deviations from sphericity 
of the clusters at liquid-state temperatures, reflecting more diversity of clusters conformations and 
shapes in order to accomodate the charges at the surface at these temperatures. 

The main effect of the addition of polarisation in the interactions between the particles in the 
MCPT calculations is to slightly decrease the number of charged particles located at the surface, 
as discussed in the previous section and shown in Fig. |2] or in other words the charged particles 
are slightly more prone to occupy the interior of the cluster. The dotted red curves in Fig. |3]a-c) 
of the Coulomb energy obtained at liquid-state temperatures of the clusters reveal that the smaller 
clusters Afl + and A\®q have intermediate behaviours between the two ideal cases, ie spherical 
conductor and sphere uniformly charged in volume, while A 7 ^ clusters depart less from the ideal 
case of the spherical conductor. 

Interestingly, the electrostatic energy curves U* MC for clusters simultaneously carrying positive 
and negative charged particles, as reported in Figure |3}i for A^jt" 14 ~, now show that these clus- 
ters behave like ideal spherical conductors for both solid-state and liquid-state temperatures. This 
certainly results from the fact that positively charged particles now, on average, solvate the neg- 
atively charged particles embedded within the volume of the clusters, thus reducing the apparent 
total charge in the interior of the cluster. The remaining positively charged particles located at the 
clusters surface then provide the final electrical thin layer of the conductor. These conclusions are 
not changed when including polarisation (red dotted curve). 
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When the cluster size is increased up to N = 1169, the amount of charged particles located 
inside the cluster increases (see Fig. |2]). For A^gg and A*?gg clusters the fraction of charged 
particles located at the surface is in the range 0.4 — 0.6, thus providing a balance between the 
number of charged particles located at the cluster surface and inside the cluster. The Coulomb 
energies U* MC for A\^(q* = 0.5) and A^g(q* = 0.2) are plotted in Fig. Hk-b and are again 
compared to the ideal Coulomb energies U* urf and U* ol . One can immediately see that the curves 
extracted from the MC calculations follow closely the case of spheres uniformly charged in volume 
over the whole range of distances r*. This behaviour is observed at both solid- and liquid-state 
temperatures. 

We thus find three general behaviours for the charged clusters investigated here: 

1. The smaller clusters (N < 147) rather behave as ideal spherical conductors at solid- 
state and liquid-state temperatures. This is consistent with the fact that the charged particles 
are predominantly located at these clusters surface for this size range. Including polarisation 
does not change these conclusions. 

2. Mixed clusters Ap~ m ~ composed of n positive and m negative charged particles behave as 
ideal spherical conductors, at both solid- and liquid-state temperatures. In such clusters, 
the negative charged particles are more prone to reside inside the clusters than the positive 
particles, while the positively charged particles make a compromise between solvating the 
negative particles within the cluster and reaching the cluster surface. As a result, these 
clusters appear as if they would carry their excess charge at the surface. 

3. The larger A^ + clusters (N = 1169) behave as ideal spheres uniformly charged within 
their volume, although the number of charged particles located at the surface roughly equals 
the number of charged particles embedded inside the cluster. Also note that locating half 
the number of charged particles at the droplets surface is not sufficient for the droplets to 
experiment some conductor-like behaviour. This is now clearly different from Konermann's 
observations 19 on a droplet composed of roughly the same number of particles (1248 in their 
study). The main difference between our two studies is the absence of a permanent dipole 
moment attached to the particles in our droplets, as well as the absence of hydrogen bonds 
between the particles. There is thus no "alignment" of dipole moments, and therefore no 
charge-dipole induced orientation of these particles in our case. Note that we conceive that 
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AiikT would favour the transfer of the excess charge onto the cluster surface, following the 
conclusions from point 2 above. 

Charged particles in our model do not carry permanent dipoles but they can carry induced 
dipoles resulting from polarisation effects. We report in Fig. 0]1 the average distribution of cos 9 
between the induced dipoles and the vector joining the cluster center of mass for A^g (<?* = 0.5) 
and Af?gg (q* = 0.2) clusters at liquid-state temperatures (See pattern in Fig. Hfc for an illustration). 
These distributions are calculated a posteriori and are averaged over 2 x 10 3 Monte Carlo con- 
figurations initially sampled without polarisation. The distributions are peaked at about cos 9 w 1 
but their tail can spread up to cos6> < (ie, 9 > 90°) for A^g (q* — 0.5) clusters. The larger 
the number n of charged particles, the sharper the distribution peaks and the shorter the extent 
of these distributions. This reflects the differences in the induced electric fields created by the 
charged particles (due to different charge distributions within the cluster) which result into more 
directionality of the induced dipoles as n increases. A finer analysis of the structure in layers of 
these clusters (not shown here) also reveals that the induced dipoles are stronger and more aligned 
along the center of mass vector for particles farther from the cluster center of mass. 

IV. CONCLUSION 

Monte Carlo simulations have been performed on A*^ + clusters composed of iV particles 
(55 < N < 1169) among which n are positively charged with charge q*, the cluster having a 
total net charge Q* = nq* . We have also investigated A^ +m ~ clusters carrying n positively and 
m negatively charged particles for a total net charge Q* = (n — m)q*. These clusters have been 
modelled with a simple and generic representation of the interactions between individual particles 
(Lennard- Jones and Electrostatic interactions, including polarisation), in dimensionless units. The 
temperature range of the simulations goes from solid-like to liquid-like states, encompassing the 
melting transition. We have investigated the localisation of the charged particles within the clus- 
ters, i.e. at the surface or embedded within the volume of the clusters, as well as the behaviour of 
these charged clusters in terms of ideal conductors or ideal spheres uniformly charged in volume. 
Both questionings stem from Electro sprayed droplets produced in Mass Spectrometry experiments 
(ESI-MS) for which these properties might be pivotal for the final production of the ions in the gas 
phase. 

For the charged clusters investigated here, with a total net positive charge Q*, the charged 
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particles tend to be predominantly located at the clusters surface at finite temperature. We found 
that an increase of Q* at fixed n or a decrease of n at fixed Q* favoured the migration of the 
positive particles to the cluster surface although complete occupancy of the surface sites was found 
out of reach for the larger A^ m clusters, even when the number of charged particles was larger 
than the number of surface sites. For the smaller clusters N < 147, low temperatures were 
found to favour the exploration of regions around the global minimum attraction basins. These 
regions are entropically unfavoured as temperature increases, especially in the melting- and liquid- 
state regimes of interest for experiments such as ESI. The localization of charged particles on the 
cluster surface still remains favoured in these regimes, with more than 70-80% occupancy for small 
clusters. Inclusion of polarisation slightly reduces the filling up of the cluster surface with charged 
particles and, interestingly, the study of charged clusters composed of both positive and negative 
particles revealed that positive particles were subject to two competing localization mechanisms: 
filling up of surface sites and solvation of negative particles. 

Considering the electrical behaviour of these charged clusters, we have found two main be- 
haviours. The smaller clusters (N < 147) behave as ideal spherical conductors at solid-state 
and liquid-state temperatures, which is consistent with the charged particles being predominantly 
located at the cluster surface. Interestingly, for the same range of clusters sizes (N = 147 here), 
mixed A^ +m ~ clusters composed of n positively and m negatively charged particles display a 
reduced number of charges located at the surface: negative charged particles are more prone to 
reside inside the clusters than the positive particles, while the positively charged particles make 
a compromise between solvating the negative particles within the cluster and reaching the cluster 
surface. Despite the reduced number of charges located at the surface of A^ +m ~, these clusters 
behave as ideal spherical conductors, where the excess charge is located within a thin layer at the 
clusters surface. On the contrary, the larger clusters (N = 1169 here) behave as ideal spheres 
uniformly charged within their volume, although the number of charged particles located at the 
surface roughly equals the number of charged particles embedded inside the cluster. 

This last result is clearly different from Konermann's results 19 on a droplet composed of roughly 
the same number of particles, i.e. 1248 microscopic water molecules (SPC model) and 10 Na + 
ions. The main difference between our two studies is the level of representation, i.e. microscopic 
representation of the molecules and atomic ions in Ref 19 and a coarse-grained representation in 
our present investigation. In our representation, there is no permanent dipole moment attached 
to the particles as well as no hydrogen bonds between the particles, contrary to the microscopic 
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representation of Konnerman and coworker. There is thus no possibility for "alignment" of dipole 
moments, and therefore no possibility for charge-dipole induced orientation of these particles in 
our droplets, which is the main argument put forward by Konermann and coworker in Ref. 19 for 
explaining the excess charge at the surface of their droplets. Note that the induced dipole moments, 
taken into account in our calculations on the A^ 9 clusters (and not considered in Ref.—), are 
not large enough to induce similar charge-dipole alignments. As a result, the only droplets 
that behave as ideal conductors are the ones composed of at most about 100 particles, because 
the charged particles are predominantly located at the clusters surface. The larger droplets of 
more than 1000 particles behave as spheres uniformly charged in volume, although roughly half 
of the charged particles are located at the clusters surface. Note that our current coarse-grained 
representation of the solvent in the droplets is very similar to the coarse-grained models of solvent 
developed and recently reviewed by Klein et al.—, i.e. chargeless and dipoless coarse grains which 
can be mapped to up to four individual microscopic water molecules. This is thus a reasonable 
model for treating coarse-grained droplets, and the results obtained here are therefore deemed 
representative of charged droplets containing coarse-grained water. 

Our present findings on charged droplets should be taken into consideration when investigat- 
ing charged droplets produced by ESI. In particular, the common assumption according which 
positive charges lie on the droplet surface should be taken with caution since the thermodynamic 
and electrical properties of the droplet may strongly affect the propensity of positive charges for 
surface localisation, as shown here. It is not clear how the electrical properties found here on 
the smaller charged droplets being ideal conductors while larger charged droplets are spheres uni- 
formly charged in volume translate for the fragmentation mechanisms of these droplets in ESI 
experiments. We are currently investigating the consequences of these findings on the fragmenta- 
tion mechanisms and possible correspondence with the Charge Residue Model (CRM) or the Ion 
Evaporation Model (IEM). We are also investigating the effects of applying an external electric 
field. 

More realistic interaction potential models should also be developped, typically inhomoge- 
neous spherical particles (with different u and a parameters), non-spherical particles, including 
patchy particles and in particular patchy proteins, and possibly going towards coarse-grained pro- 
teins embedded in droplets using the Martini or other modern coarse-grained models. This is 
another direction we are exploring. 
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VI. APPENDIX: COULOMB ENERGIES IN REDUCED UNITS 

The Coulomb energy U vo i between a probe particle with charge q pro be and the electrostatic 
potential V voi created by a sphere of radius R uniformly charged in volume with density of charge 
PQ = Q/V (Q being the total charge of the sphere and V its volume) is U vo i(r) = q pr obeV vo i(r) 
where 



In reduced units, r* = r/a for distances, U* = U/u for energies, and q* = q/y/4ire cru for 
charges, this leads to 



To properly define V* ol in reduced units, we use the simple and general expression of the electro- 
static potential created by a charge q 




r>R : V vol {r) 



pqRI 



3e r 





when posing V* ol = £ = V vol 



We then obtain 



u 




r*>R* : V v * ol = 



►*3 




3r* 
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In reduced units the Coulomb energy writes U voi = q wo beVvoi = u q^. obe V* ol = u U* ol when 
naturally posing U* ol = q* robe V* ol = U vo i/u since U voi is an energy. We eventually find 



^g; robe p* Q ( 2 



r 



*2 



r* < K ■■ t& = P ; " " -yj (6) 
r* > K : tC, = A ^ P } Rf . (7) 

For a spherical conductor with density of charge oq = Q/S (S being the surface of the sphere) 
the derivation of the Coulomb energy U* ur .j is equivalent. The expression of the electrostatic 
potential V sur f is 



r < Rq : V surf = 



t -> R ■ V - aQR ° 
r > Kq . V sur f 



e r 



which leads to 



r*<R* : U* surf = 4irq; robe a* Q R* 
_ ^q; robe a* Q Rf 



r*>R* : U* surf 
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FIG. 1: Distributions of AV*, the energy difference between the local minima explored along the 
MCPT simulations and the corresponding global minimum for selected Ag° + and A™^ clusters at 
solid-state (cyan curves), melting (black curves), and liquid-state (red curves) temperatures. For 
clarity, the breadths of Ag° + distributions have been multiplied by 2.6. 

FIG. 2: Average fraction F sur f of charged particles located at the cluster surface as a func- 
tion of the total cluster charge Q* for A^ + (Q* = nq*) clusters (open symbols) and A^ +m ~ 
(Q* = (n- m)q*) clusters (closed symbols) of size (a) N=55, (b) N=100, (c) N=147, and (d) 
N=1169. The symbols correspond to fractions obtained at the melting temperature. The error 
bars indicate the range of values spanned by F sur f for the temperatures of the Parallel Tempering 
simulations. Simulations without polarisation (black symbols), with polarisation (red symbols). 
Note that only charged particles were chosen polarisable for A™^ and A" 4 t 7 m ~ clusters. 

FIG. 3: Average Coulomb energy between a probe particle with charge \q*\ and the charged 
particles composing selected (a) A 3 5 ° 5 + , (b) A^, (c) A^, and (d) A^ 14- clusters. These 
energies are determined at solid-state temperature without polarisation (cyan solid curves) and 
liquid-state temperature without polarisation (red solid curves) and with polarisation (red dotted 
curves). Comparison with the ideal Coulomb energies U* ur f and U* ol are also reported (dashed 
curves). See text for the formula of the ideal Coulomb energies U* ur j and U* ol . 

FIG. 4: Average Coulomb energy between a probe particle with charge \q*\ and the charged par- 
ticles composing selected (a) A^gg and (b) Af^g clusters. Comparison with the ideal Coulomb 
energies U* ur f and U* ol are also reported (dashed curves), (c) Pattern that defines the angle 9 used 
for plotting cosine distributions P CO s6» of large clusters, (d) Average P cos e distributions for the 
estimated induced dipoles on Aj™g (q* = 0.5) (black curves) and A 4 ?^ (q* = 0.2) (red curves). 
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